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Abstract.  Problems  involving  the  management  of  groundwater  resources  occur  routinely, 
and  management  decisions  based  upon  optimization  approaches  offer  the  potential  to  save 
substantial  amounts  of  money.  However,  this  class  of  application  is  notoriously  difficult  to 
solve  due  to  non-convex  objective  functions  with  multiple  local  minima  and  both  nonlinear 
models  and  nonlinear  constraints.  We  solve  a  subset  of  community  test  problems  from  this  ap¬ 
plication  field  using  MODFLOW,  a  standard  groundwater  flow  model,  and  IFFCO,  an  implicit 
filtering  algorithm  that  was  designed  to  solve  problems  similar  to  those  of  focus  in  this  work. 
While  sampling  methods  have  received  only  scant  attention  in  the  groundwater  optimization 
literature,  we  show  encouraging  results  that  suggest  they  are  deserving  of  more  widespread 
consideration  for  this  class  of  problems.  In  keeping  with  our  objectives  for  the  community 
problems,  we  have  packaged  the  approaches  used  in  this  work  to  facilitate  additional  work 
on  these  problems  by  others  and  the  application  of  implicit  filtering  to  other  problems  in  this 
field.  We  provide  the  data  for  our  formulation  and  solution  on  the  web. 

Keywords:  Implicit  filtering,  Well  field  design.  Groundwater  flow  and  transport 


1.  Introduction 

Groundwater  resources  are  important  because  about  50%  of  the  population 
of  the  United  States  rely  upon  this  resource  for  drinking  water.  Typical  goals 
in  groundwater  resources  management  include  designing  and  managing  sys¬ 
tems  to  supply  drinking  water  and  to  restore  contaminated  drinking  water  to 
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a  usable  quality  at  a  minimum  cost.  Accomplishing  these  goals  requires  a 
model  to  describe  the  system  of  concern,  an  appropriate  objective  function, 
constraints,  and  an  optimization  algorithm.  The  objective  function  and  con¬ 
straints  provide  the  linkage  between  the  simulation  model  and  the  optimizer. 
Because  porous  media  systems  arc  typically  heterogeneous  over  small  scales 
and  described  by  nonlinear-  processes,  subsurface  simulators  can  be  expensive 
to  evaluate  and  subject  to  uncertainty,  or  stochastic  in  nature.  The  resulting 
optimization  problems  can  also  be  diffi  cult,  with  objective  functions  that  are 
non-convex  and  have  multiple  local  minima,  and  both  models  and  constraints 
that  are  nonlinear.  For  these  reasons,  subsurface  optimization  problems  are 
both  important  and  challenging. 

To  aid  the  evolution  of  optimal  design  of  subsurface  flaw  and  transport 
applications,  a  set  of  community  problems  (CP’s)  was  developed,  [30],  that 
are  typical  of  problems  commonly  encountered  and  which  cover  a  range  of 
complexity.  It  was  reasoned  that  focusing  on  a  common  set  of  CP’s  would 
allow  for  not  only  advancement  of  approaches  to  solve  an  important  set  of 
problems,  but  also  a  means  to  aid  comparison  of  various  aspects  of  the  so¬ 
lution  approach  on  the  same  set  of  problems.  It  was  also  anticipated  that  the 
CP’s  would  catalyze  the  introduction  of  new  classes  of  optimization  methods 
into  the  groundwater  fi  eld  and  result  in  more  active  participation  of  the  ap¬ 
plied  mathematics  community  in  the  evolution  of  solution  approaches  for  this 
class  of  application. 

A  subset  of  the  CP’s  is  a  standard  water  supply  application.  Roughly 
speaking,  the  objective  is  to  locate  a  set  of  water  supply  production  wells  and 
fi  nd  their  pumping  rates  such  that  cost  is  minimized  subject  to  constraints  on 
the  total  amount  of  water  that  must  be  produced,  the  hydraulic  head  in  the 
wells,  the  production  capacity  of  a  well,  and  the  portion  of  the  domain  over 
which  a  well  may  be  located.  Evaluation  of  the  objective  function  requires 
a  groundwater  flaw  simulator  that  solves  for  hydraulic  head  as  a  function 
of  space  and  time  given  a  spatial  and  temporal  domain,  material  properties, 
auxiliary  conditions,  and  well  design  information. 

Common  approaches  for  solving  water  supply  management  problems  in¬ 
clude  gradient-based  methods  or  genetic  algorithm  (GA)  approaches  in  which 
a  candidate  set  of  well  locations  is  selected  and  the  pumping  rate  of  each  well 
is  a  design  variable.  Gradient-based  methods  are  not  reliable  because  of  the 
non-convex,  noisy  nature  of  the  problem.  Global  optimization  methods  such 
as  genetic  and  evolutionary  algorithms  [38,  40,  24,  1],  simulated  annealing 
[28,  39],  and  tabu  search  [45]  have  been  applied  to  many  subsurface  remedi¬ 
ation  problems  (see  [30]  for  many  more  references).  However,  these  methods 
can  be  expensive.  Sampling  methods  are  a  potentially  attractive  alternative 
class  of  approach  for  this  sort  of  problem,  which  have  received  scant  attention 
in  the  water  resources  community  but  have  been  concluded  to  be  deserving 
of  a  more  thorough  investigation. 
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In  this  paper  we  take  the  view  that  there  is  enough  structure  in  the  problem 
to  use  a  deterministic  sampling  method.  These  methods  arc  designed  to  solve 
problems  with  diffi  cult,  but  not  violently  oscillatory  optimization  landscapes, 
such  as  the  ones  in  Figures  8,  9,  10,  11  in  §  8.  Gradient-based  methods  arc 
likely  to  have  trouble  with  such  problems,  either  fi  nding  local  minima,  stag¬ 
nating,  or  failing  to  fi  nd  a  descent  direction.  In  our  testing  of  a  gradient-based 
method,  which  we  report  in  §  8,  we  observed  this  type  of  failure.  The  Nelder- 
Mead  [34],  Hooke-Jeeves  [23],  MDS  [16,  43],  DIRECT  [25],  and  implicit 
fi  ltering  [21,  20,  26]  are  examples  of  discrete  sampling  methods. 

The  objectives  of  this  work  are:  (1)  to  provide  an  initial  analysis  of  a  subset 
of  groundwater  CP’s,  which  have  been  recently  published;  (2)  to  formulate 
a  solution  to  these  problems  with  a  sampling  method,  in  this  case  implicit 
fi  ltering;  (3)  to  compare  the  results  with  a  genetic  algorithm  approach,  and 
to  explain  why  traditional  gradient-based  methods  can  and  do  fail;  (4)  to 
examine  the  characteristics  of  the  solution  space  and  illustrate  the  challenges 
that  this  class  of  problem  poses;  and  (5)  to  point  the  way  toward  future 
improvements  for  the  solution  of  this  class  of  problem. 


2.  Conceptual  Model 

An  aquifer  is  a  fully  saturated,  water-hearing  region  and  is  considered  con- 
fi  ned  if  bounded  on  both  the  top  and  the  bottom  by  essentially  impermeable 
material.  An  unconfi  ned  aquifer  has  the  water  table  as  its  upper  bound.  The 
main  difference  between  the  two  geological  formations  is  that  the  saturated 
thickness  of  an  unconfi  ned  aquifer  varies  as  the  hydraulic  head  varies,  thus 
leading  to  a  nonlinear  tree-boundary  problem. 

We  consider  a  well-fi  eld  design  problem.  The  hydrological  settings  arc 
homogeneous  confi  ned  and  unconfi  ned  aquifers  in  three  spatial  dimensions. 
For  the  problems  considered,  a  set  of  wells  is  distributed  in  the  domain.  Each 
well  is  allowed  either  to  inject  or  extract  water.  Well-fi  eld  design  problems 
involve  the  selection  of  well  locations  and  pumping  rates  to  minimize  the 
cost  of  water  production.  The  cost  of  supplying  water  typically  involves  the 
cost  to  drill,  equip,  and  connect  wells  to  a  treatment  or  distribution  system, 
and  the  cost  to  pump  the  water  and  maintain  the  well.  In  turn,  the  cost  to 
pump  groundwater  depends  upon  the  energy  needed  to  lift  the  water  from  its 
level  below  the  ground  surface  to  the  discharge  point  and  to  supply  suffi  dent 
discharge  pressure  to  achieve  the  desired  fbw. 

The  decision  variables  for  this  type  of  problem  are  the  pumping  rates 
{Qi}?=  i  ( m3/s )  at  the  n  wells  in  the  model  and  the  locations  {(aq,  2/i)}f=1 
of  the  wells.  Pumping  rates  can  be  constant  or  variable  in  time  depending 
upon  the  application.  In  the  application  considered  in  this  paper,  a  constant 
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few  rate  is  realistic.  This  is  because  any  transients  decay  very  early  in  the 
fi  ve  year  time  horizon. 


3.  Formulation 


The  physical  domain,  see  Figure  1,  is  17  =  [0, 1000]  x  [0, 1000]  x  [0,  30]  m 
with  the  ground  elevation  at  zgs  =  60  m  for  the  confi  ned  aquifer  and  ZgS  =  30 
m  for  the  unconfi  ned  aquifer. 

Flow  in  saturated  porous  media  can  be  described,  [30],  by 

f)h 

Ss-  =  V-(KVh)+S,  (1) 

where  Ss  (1/m)  is  the  specifi  c  storage  coeffi  cient,  the  unknown  h  (m)  is  the 
hydraulic  head,  K  (m/s)  is  the  hydraulic  conductivity  [30].  Here  the  source 
term  S  is  a  model  of  the  wells,  a  sum  of  5-functions  that  satisfi  es 


in 


S(t)d£l  =  ^2  Qi • 


i=l 


(2) 


17  is  the  spatial  domain.  The  wells  are  assumed  to  extract  from  near  the  bottom 
of  the  aquifer.  If  a  numerical  solution  is  discrete  in  the  ^-dimension  then  only 
the  bottom  layer  of  cells/elements  should  include  the  well  source  terms. 


tions: 


dh 

dx 


x=0 


dh 

dy 


y= o 


:  following  boundary 

and  initial  condi- 

dh 

(3) 

dz 

=  0,/  >  0 

2=0 

-1.903  x  10 ~8 (m/s) 

(4) 

=  50 

—  O.OOly(m) 

(5) 

=  50 

—  O.OOlx'(m) 

(6) 

0)  = 

hs 

(V) 

Here 


Tfdh 
qz  =  -K— 

dz 


is  the  Darcy  fbx  out  of  the  domain,  a  negative  sign  in  eqn  (4)  thus  represents 
fbw  into  the  aquifer  or  recharge  that  could  be  the  result  of  rainfall  inti  ltration 
or  leakage  from  an  overlying  aquifer,  and  hs  is  the  steady  state  solution  to  the 
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few  problem  prior  to  the  addition  of  wells.  We  use  S  s  =  10  6  (1/m).  For  the 
unconfi  ned  aquifer,  (4),  (5)  and  (6)  are  replaced  with 

qz(x,y,h,t  >  0)  =  —1.903  x  10~8(m/s),  (8) 

/t(1000,  y,z,t>  0)  =  20  -  O.OOly(m),  (9) 

and 

h{x,  1000,  z,t>  0)  =  20  -  O.OOlx(m).  (10) 

Ss  =  2.0  x  HP  1  is  the  specific  yield  of  the  unconfined  aquifer.  For  the 
homogeneous  applications,  K  =  5.01  x  10-5  (m/s). 
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y=  1,000  (m) 


Dirichlet  Boundary  Conditions 


o 

U 


■d 

a 

3 

O 

PQ 

* 

o 

E 

0 

Z 


Recharge  into  Aquifer  Top,z=30  - 


x=0  (m 
y=0  (m 


No  Flow  Boundary  Conditions 


h=19  or  49  (m) 
x=  1,000  (m) 
y= 1,000  (m) 


O 

£ 

cT 

Od 

o 

c 

3 

a 

n 

0 

3 

a 

o' 


h=20  or  50  (m) 
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Figure  1.  Physical  Domain 


4.  Objective  Function 

We  consider  a  capital  cost  fc  and  an  operational  cost  f°  seeking  to  minimize 
fT  =  fc  +  f°.  The  objective  function  depends  on  the  pumping  rates  {Qi}f=1 
and  locations  {(x^,  yi)}f=1  of  n  operating  wells.  Note  that  Qi  <  0  means  the 
well  is  extracting  water,  and  Qi  >  0  means  the  well  is  injecting  water.  For 
this  work,  we  begin  with  a  virtual  fi  xed  well  fi  eld  containing  N  wells  with 
the  number  of  operating  wells  n  <  N. 
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Since  there  is  a  fi  xed  installation  cost  for  wells,  an  important  aspect  of  the 
optimization  procedure  is  the  manner  in  which  wells  arc  removed  from  the 
design,  thereby  signifi  candy  decreasing  the  total  cost.  A  well  is  considered 
installed  and  operating  if  \Qi\  >  0.0001.  If  the  optimizer  specifies  a  value 
with  | Q,  |  <  0.0001  then  we  neither  apply  the  well  source  term  nor  do  we 
include  the  cost  of  the  well  in  the  objective  function.  This  approach  results  in 
non-smoothness  in  the  objective  function,  but  provides  a  reasonable  mecha¬ 
nism  for  removing  wells  from  the  design  that  our  optimizer  was  capable  of 
triggering. 

The  objective  function  is  given  by 

n 

/T  =  E^°  +  53  Cl\QT\bl(zgs-hmin)b2+  (11) 


where  the  cost  coeffi  dents  q  and  exponents  bj  arc  given  in  Table  I.  Here  d, 
is  the  depth  of  well  i,  Q™  is  the  design  pumping  rate,  hmtn  is  the  minimum 
allowable  head,  hi  is  the  hydraulic  head  in  well  i,  t  j  is  the  total  design  time, 
which  was  taken  as  5  years,  and  zgs  is  the  elevation  of  the  ground  surface. 
Injection  wells  are  assumed  to  operate  under  gravity  feed.  In  /c,  the  fi  rst  term 
denotes  the  cost  to  install  all  the  wells,  and  the  second  term  accounts  for  the 
additional  cost  for  pumps  for  the  extraction  wells.  In  f°  we  have  a  lift  cost 
that  applies  to  the  extraction  wells  and  an  injection  cost  that  applies  to  the 
injection  wells. 

The  design  pumping  rates  {Q]n},  i.e.  the  maximum  rates  at  which  a  given 
well  can  pump,  depend  upon  the  aquifer  properties,  casing  and  discharge 
piping  size,  pump  characteristics,  screen  length  and  opening  size,  effective¬ 
ness  of  the  development,  and  local  geochemical  conditions.  One  could,  in 
principal,  treat  the  properties  of  the  wells  and  pumps  as  optimization  param¬ 
eters.  We  do  not  do  this  here,  and  focus  on  more  fundamental  aspects  of  the 
formulation. 


1=1 


Q;<  — .0001 


E  c2Qi(hi  -  zgs )  +  E  C3Q1 


.0001 


j,Qi>.0001 


f° 


5.  Constraints 

We  constrain  the  hydraulic  head  and  pumping  rates  for  the  objective  function 
given  in  (1 1).  The  constraints  arc  given  by 
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Table  I.  Objective  Function  Data 


data 

value 

units 

co 

5.5  x  103 

%/mb 

Cl 

5.75  x  103 

$/[(m3/s)61  •  mb2] 

C2 

2.90  x  10-4 

$/m4 

C3 

1.45  x  10-4 

CO 

,s_ 

bo 

0.3 

- 

bi 

0.45 

- 

&2 

0.64 

- 

Zgs 

60  confined 

m 

ZgS 

30  unconfined 

m 

di 

Zgs 

m 

QT 

1.5  Qi 

m3 /s 

and 


n 


Qt  =  ^2Qi<  Qt™, 

(12) 

i= 1 

<Qi<Qimax ,  i  =  l, n, 

(13) 

<  hi  <  hmax,  i  =  1, ...,  n, 

(14) 

where  Qr  is  the  net  pumping  rate,  is  minimum  allowable  total  ex¬ 

traction  rate,  Qemax  is  the  maximum  extraction  rate  at  any  well,  Qimax  is 
the  maximum  injection  rate  at  any  well,  hrnnx  is  the  maximum  allowable 
head,  and  limin  is  the  minimum  allowable  head.  Values  for  the  bounds  in  the 
constraints  arc  given  in  Table  II.  We  require  that  the  wells  be  at  least  200  m 
from  the  boundary  on  which  Dirichlet  boundary  conditions  are  applied,  i.  e. 


0  <  Xi,  Hi  <  800. 


(15) 


In  addition  to  (15),  we  do  not  allow  two  wells  to  occupy  the  same  grid 
point.  In  the  course  of  the  optimization,  if  two  wells  converge  to  the  same 
location,  our  choice  of  simulator  would  implement  the  two  wells  as  one  well, 
operating  at  the  sum  of  the  two  pumping  rates.  In  turn,  only  one  well  would 
operate  in  the  fbw  simulation,  yet  two  wells  would  be  included  in  the  instal¬ 
lation  cost.  For  our  choice  of  spatial  discretization,  this  indirectly  implies  that 
the  distance  between  wells  is  at  least  20m  apart. 

Constraint  (12)  sets  a  minimum  target  for  extraction,  which  is  the  purpose 
of  the  well  fi  eld.  For  the  problem  considered  here,  the  installation  costs  arc  far 
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more  that  the  operating  costs  for  a  single  year.  Therefore,  once  the  minimum 
extraction  target  is  reached,  it  would  only  make  sense  to  drill  additional  wells 
if  the  long-term  operating  savings  is  signifi  cant.  Since  fi  ve  wells  extracting  at 
the  maximum  level  satisfy  (12)  with  equality,  one  logical  formulation  of  the 
problem  is  to  fi  nd  the  optimal  location  of  fi  ve  wells,  each  extracting  as  much 
as  possible. 

Constraint  (13)  refbcts  physical  limits  on  the  pumps  and  well  design.  Well 
designs  arc  typically  limited  by  the  size  distribution  of  the  porous  medium 
and  the  resulting  size  of  the  well  screen. 

The  upper  bound  in  constraint  (14)  keeps  the  hydraulic  head  below  the 
surface  elevation  and  the  lower  bound  ensures  that  excessive  drawdown  will 
not  occur.  This  constraint  is  a  linear  function  of  the  pumping  rates  for  the 
confi  ned  case  but  a  nonlinear  function  for  the  unconfi ned  case,  and  in  both 
cases  a  highly  nonlinear  function  of  the  locations  of  the  wells. 

Table  II.  Constraint  Data 


data 

value 

units 

z- \min 

-3.2  x  10“2 

m3 /s 

Qemax 

-6.4  x  10“3 

m3 /s 

Qimax 

6.4  x  10~3 

m3 /s 

j^rnin 

40  confined 

m 

j^rnax 

60  confined 

m 

yrnin 

10  unconfined 

m 

j^rnax 

30  unconfined 

m 

5.1.  Optimization  Problem  Formulation 

In  this  section  we  describe  how  we  packaged  the  problem  for  the  optimization 
algorithm.  The  objective  function  f1  is  discontinuous,  and  some  of  the  con¬ 
straints  (13)  and  (15)  are  simple  bounds  on  the  variables.  Implicit  filtering, 
the  optimization  method  we  use  in  this  paper,  is  designed  to  handle  diffi  cult 
objective  functions  and  bound  constraints. 

If  we  set  n  =  5,  then  the  constraints  (12)  and  (13)  require  the  pumping 
rates  to  be  exactly  Q™m/ 5.  Thus,  in  this  situation,  we  need  only  optimize 
well  locations  and  apply  constraint  (14).  If  we  set  n  >  5,  then  we  must  also 
optimize  pumping  rates  and,  therefore,  all  the  constraints  must  be  enforced  by 
the  optimizer.  Constraint  (14)  is  highly  nonlinear  while  constraint  (12)  is  not  a 
box  constraint,  and  neither  constraint  can  be  handled  directly  by  the  projected 
quasi-Newton  algorithm.  In  this  case  the  objective  function  returns  a  failure 


katie.tex;  11/07/2003;  8:20;  p.8 


9 


when  either  (12)  or  (14)  arc  violated.  Our  implementation  of  implicit  fi  ltering 
will  assign  an  artifi  cial  (see  §  6.2)  value  to  the  function  when  it  returns  a 
failure.  This  is  a  standard  approach  for  handling  nonlinear  constraints  in  many 
sampling  methods,  [42,  9,  27]. 

We  will  fi  x  the  number  of  wells  and  consider  the  vector  of  design  variables 

Z  —  (xi  Xn ,  2/1 ,  ...  ,  2/nj  Qli  ■  ■  ■  i  Qn)  £  R 


We  defi  ne  the  feasible  set  for  the  bound  constraints  as 

V0  =  {Z\  (13)  and  (15)  hold.}  =  {Z  \  Z™in  <  Z>  <  Z™ax}.  (16) 
Our  optimization  problem  is 

S/T(z>-  (17> 

where  f1  is  given  by  (1 1)  if  (12)  and  (14)  arc  satisfi  ed  and  a  failure  is  reported 
if  either  of  (12)  or  (14)  arc  violated. 


6.  Implicit  Filtering 

The  objective  function  is  highly  nonlinear  and  non-convex,  discontinuous 
because  of  the  jumps  as  wells  arc  added  and  deleted,  and  noisy,  because  of 
internal  iterations  in  the  simulators.  For  these  reasons,  as  we  said  in  §  1, 
a  conventional  gradient-based  optimization  method  may  fail.  A  sampling 
method,  which  only  evaluates  the  objecive  function  and  constraints  to  guide 
the  optimization,  is  most  appropriate  for  this  kind  of  problem. 

In  this  paper  we  use  IFFCO  [10],  a  FORTRAN  implementation  of  the 
implicit  filtering  algorithm  [26,  21,  20].  We  based  this  decision  on  our  own 
familiarity  with  the  optimizer  and  our  past  success  with  it  on  other  problems 
of  a  similar  mathematical  nature  [4,  42,  9],  although  we  arc  not  aware  of 
any  use  of  implicit  fi  ltering  for  the  type  of  application  problem  of  concern 
in  this  work.  This  choice  signifi  candy  in  flic  need  the  decisions  on  handling 
constraints  and  the  locations  of  the  wells. 

Implicit  fi  ltering  has  been  described  in  detail  and  analyzed  elsewhere.  We 
refer  the  reader  to  [26]  for  the  details  of  the  algorithm  and  to  [26,  21,  11]  for 
convergence  analysis.  In  §  6.1  we  sketch  the  algorithm  and  its  implementa¬ 
tion  in  IFFCO  only  in  enough  detail  to  explain  how  this  choice  affected  the 
formulation  of  the  problem. 
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6.1.  The  ALGORITHM 

Implicit  fi  ltering  is  a  projected  quasi-Newton  method  that  uses  fi  nite  dif¬ 
ference  gradients.  The  difference  increment  is  reduced  as  the  optimization 
progresses,  thereby  avoiding  some  local  minima,  discontinuities,  or  nons¬ 
mooth  regions  that  would  trap  a  conventional  gradient-based  method.  The 
problems  considered  in  this  paper  arc  exactly  the  kind  that  the  method  was 
designed  to  solve. 

Implicit  filtering  begins  by  rescaling  the  variables  so  that  the  feasible 
region  is 

£>  =  {£|0  <  &  <  1}.  (18) 

We  will  discuss  the  algorithm  in  terms  of  the  scaled  feasible  region  in  (18) 
but  the  application  in  terms  of  the  actual  bounds  (16). 

To  make  the  transition  from  fT  to  the  scaled  form,  we  defi  ne  £  component¬ 
wise  by 

&  =  (Zi  -  z™in)/(Z™ax  -  Zfn) 

and  let 

no  =  fT(z). 

The  optimization  problem  for  /  is  now 

min^evf(0- 

For  a  given  difference  increment  (called  a  scale)  S  G  (0, 1/2]  and  (eP, 
we  let  VsfiO  be  the  difference  gradient  whose  components  arc 

—  the  central  difference  gradient  in  the  zth  coordinate  direction  if  both  of 

£  ±  Sei  G  V,  or 

—  the  one-sided  difference  gradient  in  the  i  coordinate  direction  if  only  one 

of  £  ±  Sei  G  £>■ 

Since  S  <  1/2,  at  least  one  of  £  ±  Sei  G  £>.  We  let  the  stencil  ,S'(£)  be 
those  points  in  the  centered  difference  stencil  that  arc  in  V  and  used  in  the 
computation  of  Vsf-  If 

/(£)  <  min  f(rf)  (19) 

ves  (0 

we  say  that  stencil  failure  has  occurred  and  terminate  the  quasi-Newton 
iteration  at  that  scale. 

If  H  is  a  model  Hessian,  a  projected  quasi-Newton  iteration  from  £  has 
the  general  form 

£(A)=)P(£-A 
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where  V  is  the  projection  onto  V 

(0  if  Ci  <  o 

V(Oi  =  \  Ci  if  0  <  ^  <  1 

U  if  e*  >  i 

In  IFFCO,  the  step  length  A  is  computed  with  a  quadratic  model  [10]  and 
a  step  is  accepted  if  the  suffi  cient  decrease  condition 

m  -  m  <  a\75f(z)Tm)  -  o,  (20) 

holds.  In  IFFCO,  as  is  standard,  a  =  10-4.  We  say  that  the  quasi-Newton 
iteration  is  successful  if 

||£-£(1)||<t<5.  (21) 

The  algorithmic  parameter  r  can  have  a  signifi  cant  effect  on  the  performance 
of  the  optimization.  For  the  problems  we  consider  here,  however,  we  were 
able  to  successfully  use  the  default  value  of  r  =  1. 

The  fi  nite  difference  projected  quasi-Newton  loop  in  IFFCO  is  summa¬ 
rized  in  algorithm  fdquasi.  fdquasi  is  a  naturally  parallel  algorithm;  all 
the  function  evaluations  needed  to  compute  V sf  can  be  done  in  parallel.  We 
exploited  this  simple  parallelism  to  perform  the  computations  reported  in  this 
paper. 


Algorithm  1  fdquasi  (£,  f.  prnax,  r,  S,  amax) 

p  =  1 

while  p  <  pmax  and  ||£  —  V(£  —  V,5/(£))||  >  tS  do 
compute  /  and  Vsf 
if  (19)  holds  then 

terminate  and  report  stencil  failure 
end  if 

update  the  model  Hessian  H  if  appropriate;  solve  Hd  =  —VsfiO 
use  a  backtracking  line  search,  with  at  most  amax  backtracks,  to  fi  nd  a 
step  length  A 

if  amax  backtracks  have  been  taken  then 
terminate  and  report  line  search  failure 
end  if 

a^v(t+\  d) 

p  p + 1 

end  while 

if  p  >  pmax  report  iteration  count  failure 


Implicit  fi  ltering  calls  fdquasi  repeatedly  with  a  sequence  of  scales  {41- 
Algorithm  imfilter  is  a  simple  sketch. 
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Algorithm  2  imfilter(£.  /,  pmax,  r,  {^},  amax) 

for  k  =  0, ...  do 

f  dquasi(£,  f,pmax,  r,  amax) 

end  for 


The  algorithmic  parameters  that  arc  important  to  implicit  fi  ltering  arc  the 
limit  amax  on  the  number  of  step  size  reductions,  pmax  on  the  number  of 
nonlinear  iterations,  and  the  parameter  r  in  the  termination  criterion.  For  the 
calculations  reported  here,  we  set  pmax  =  100  (the  default),  r  =  1  (the 
default),  and  amax  =  3  (the  default).  The  parameters  in  imfilter  that  control 
the  quasi-Newton  loop  arc  the  sequence  of  scales  {/>/-}■  Our  choice  in  this 
work  was 

Sk  =  2~k-\  0  <  k  <  10. 

The  analysis  of  implicit  fi  ltering  begins  with  the  paradigm 

f  =  fs  +  (/>  (22) 

where  fs  is  a  smooth  function  and  <i>  represents  the  “noise”  in  the  problem. 
For  the  theoretical  convergence  results  in  [42,  26,  1 1]  we  assume  that  6  is  an 
everywhere-defi  ned  function  on  O  and  set 

=  max  l<M£)l- 

f?eS(0 

One  can  show  that  if  either  (19)  or  (21)  hold,  that 

r(^-V/s(0)ll=O(«J+Ms(0/<5).  (23) 

The  convergence  theory  for  implicit  fi  ltering  [26,  21,  11]  arc  based  on  (23). 

IFFCO  supports  the  SRI  [7,  18]  and  the  BFGS  [41,  8,  19,  22]  quasi- 
Newton  models  of  the  Hessian.  We  used  the  SRI  update  in  this  paper.  In  our 
experience  the  SRI  update  performs  better  for  bound-constrained  problems. 

Implicit  fi  ltering  can  be  restarted  after  it  terminates  and  the  convergence 
theory  [21]  is  stronger  if  one  does  that.  In  practice,  restarting  usually  has  no 
effect.  For  the  problems  in  this  paper,  however,  we  had  to  restart  IFFCO  once 
to  obtain  consistently  good  results. 

6.2.  Failure  of  the  function 

IFFCO  responds  to  a  failure  of  /  in  two  ways.  If  the  failed  function  evaluation 
fT(z )  is  part  of  the  evaluation  of  V  hfT(Q>  then  an  art  ill  cial  value  [9]  of 

r  +  io-6in 
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is  assigned  to  f(z).  Here  /*  is  the  largest  function  value  in  the  stencil  S(£).  If 
the  function  evaluation  failure  is  part  of  the  line  search,  the  the  value  f  scale 
is  assigned  to  fT. 

f  scale  is  an  approximation  to  the  maximum  value  of  / 7  in  the  feasible 
set  Vq  for  the  bound  constraints  (16).  We  set  f  scale  to  20%  more  than  the 
value  of  /  at  the  initial  iterate  in  this  paper. 

This  approach  to  handling  constraints  is  natural  if  the  failure  of  the  objec¬ 
tive  function  is  a  consequence  of,  for  example,  an  internal  iteration’s  failure  to 
converge.  In  the  case  of  the  problem  considered  here,  while  the  constraints  arc 
directly  specili  ed  by  (14),  the  evaluation  of  /;,  requires  a  call  to  the  simulator 
which,  as  a  function  of  the  well  locations,  is  highly  nonlinear  even  for  the 
continuous  problem.  For  the  discrete  problem  considered  here,  where  the 
well  locations  arc  rounded  to  grid  points  before  the  call  to  the  simulator,  the 
constraint  function  is  discontinuous. 


7.  Evaluation  of  the  Objective  Function 


IFFCO  requires  an  external  subroutine  to  evaluate  the  objective  function  fT. 

To  do  this  we  must  compute  the  hydraulic  head  values,  {hi},  at  the  well 
locations  {(xi,yi)}  for  a  given  set  of  pumping  rates  {Q,}-  Computation  of 
{hi}  uses  a  groundwater  fbw  simulator  to  solve  (1).  For  this  work  we  use  the 
U.S.  Geological  Survey  code  MODFLOW-96  [31].  MODFLOW  is  a  block- 
centered  fi  nite  difference  code  that  simulates  saturated  groundwater  fbw  and 
allows  for  a  variety  of  boundary  conditions  and  irregular  physical  domains. 
MODFLOW  is  widely  used  and  well  supported. 

A  MODFLOW  simulation  requires  an  input  fi  le  containing  the  location 
and  pumping  rates  of  the  wells  in  the  model.  If  n  >  5,  each  function  evalu¬ 
ation  requires  a  new  set  of  pumping  rates  and  thus  the  MODFLOW  well  fi  le 
must  be  created  each  time  the  objective  function  is  evaluated.  Moreover,  once 
the  MODFLOW  simulation  is  complete,  the  values  of  h,  must  be  extracted 
from  the  MODFLOW  output  fi  le.  A  typical  function  evaluation  is  shown  in 
Figure  2. 

To  generate  the  necessary  data  fi  les  to  run  MODFLOW  we  used  the  Groundwater 
Modeling  System  (GMS),  version  3.1.  GMS  is  a  modular  interface  to  a  vari¬ 
ety  of  fbw  and  transport  codes,  including  MODFLOW.  GMS  has  a  graphical 
environment  that  allows  the  user  to  generate  grids,  defi  ne  characteristics  of 
the  porous  media,  and  visualize  solutions.  GMS  was  used  to  generate  the 
starting  heads  for  (7),  to  create  the  necessary  data  fi  les  for  MODFLOW,  to 
determine  an  appropriate  initial  iterate  for  the  optimization,  and  then  again  to 
test  the  results  of  the  optimizer. 
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Round  Well  Locations 
To  Grid  Points 


Write  Pumping  Rates 
and  Well  Locations 
To  MODFLOW  Well  File 


Run  MODFLOW 
Simulation 


Figure  2.  Objective  Function  Evaluation 


8.  Numerical  Results 

We  consider  two  formulations  of  the  design  space.  For  this  work,  the  instal¬ 
lation  cost  of  an  extraction  well,  which  is  roughly  $20,000,  is  high  compared 
to  the  annual  operating  cost  which  is  roughly  $1,000.  Since  n  =  5  wells 
extracting  at  Qemax  =  — 0.0064(m3/s)  satisfies  the  water  supply  demand 
(12)  exactly,  one  obvious  formulation  is  to  fi  x  n  =  5  and  {Q*}f=1  =  -0.0064 
and  seek  the  optimal  locations,  {(xj,  r/*)}f=1  to  minimize  only  the  operational 
cost  (f°  in  (11)).  We  also  include  a  formulation  in  which  we  start  with  N  =  6 
wells  and  seek  the  optimal  locations  and  pumping  rates  to  minimize  (11), 
fT  =  fc  +  f°.  Intuitively,  since  the  installation  cost  is  so  high,  we  would 
expect  the  fi  ve  well  confi  guration  to  have  the  lowest  cost. 

We  examined  the  performance  of  one  gradient-based  code,  the  FDNIPS 
solver  from  the  OPT++  v2.0  [33]  framework.  This  code  is  a  nonlinear-  interior 
point  code  based  on  the  work  in  [17,  2,  3].  The  code  uses  finite  difference 
gradients,  either  trust  region  or  line  search  globalization,  and  a  choice  of  three 
merit  functions.  We  tried  several  combinations  of  the  options.  In  every  case 
the  optimization  failed  after  1000  calls  to  the  function  or  failed  because  the 
line  search  had  reduced  the  step  length  40  times  without  a  suffi  cient  decrease 
in  the  merit  function. 

For  comparison,  we  include  results  obtained  with  a  simple  genetic  algo¬ 
rithm  (GA).  The  performance  of  a  GA  for  constrained  optimization  problems 
often  depends  strongly  on  a  number  of  factors  including  the  method  used 
to  encode  the  design  variables,  the  choice  of  selection,  crossover,  and  mu- 
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tation  operators,  and  the  manner  in  which  constraints  arc  handled  [36,  29]. 
Here,  we  consider  a  single-objective  GA  which  incorporates  both  real-  and 
binary-coded  variables,  and  uses  binary  tournament  selection  [15],  For  the 
real-coded  variables,  the  simulated  binary  crossover  (SBX)  operator  [15,  14] 
with  polynomial  mutation  is  used  while  single -point  crossover  with  bitwise 
mutation  are  used  for  binary-coded  variables. 

An  approach  based  on  [12]  is  used  to  include  constraints  without  the  use 
of  penalty  parameters.  Box  constraints  such  as  those  in  (13)  arc  enforced 
automatically  in  the  generation  of  candidate  design  variables,  while  a  con¬ 
straint  such  as  (12)  is  formulated  as  a  non-negative  function  g(Z)  >  0.  The 
GA  tournament  selection  process  is  then  modifi  ed  to  account  for  the  three 
scenarios:  (1)  when  two  feasible  solutions  arc  compared,  the  one  with  lower 
objective  value  is  preferred;  (2)  when  a  feasible  and  infeasible  solution  arc 
compared,  the  feasible  one  is  taken;  and  (3)  when  two  infeasible  solutions  arc 
compared  the  one  with  lower  overall  constraint  violation  is  preferred  [12], 

Parameters  like  the  population  size,  number  of  generations,  as  well  as  the 
probabilities  and  distribution  indexes  chosen  for  the  crossover  and  mutation 
operators  effect  the  performance  of  a  GA  [36,  29],  For  the  purposes  of  our 
comparison,  we  wished  to  limit  the  number  of  simulations  performed  by 
the  GA  to  a  range  of  2  to  3  times  the  number  required  by  IFFCO.  Since 
the  total  number  of  objective  function  evaluations  is  roughly  the  product  of 
the  population  size  and  number  of  generations,  this  restricted  our  choices  to 
fairly  small  populations  and  few  generations.  We  also  wished  to  use  similar 
parameter  values  across  the  various  problems.  Although  we  did  not  perform 
a  systematic  study  to  fi  nd  the  best  possible  combinations,  we  experimented 
with  a  series  of  population  sizes,  numbers  of  generations,  and  crossover  and 
mutation  parameters  to  fi  nd  a  combination  that  gave  representative  perfor¬ 
mance  for  each  of  the  problems  we  considered.  Unless  noted,  the  values  used 
are  listed  Table  III. 

The  GA  code  used  is  implemented  in  C  and  is  available  for  download 
from  [13].  The  user  is  required  to  implement  problem-specific  routines  for 
evaluating  objective  functions  and  constraints. 
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Table  III.  GA  parameters 


30  size  of  population 

30  number  of  generations 

0.9  crossover  probability 

0. 1  real-coded  mutation  probability 

20  distribution  index  for  real-coded  crossover 

10  distribution  index  for  real-coded  mutation 

0.5  binary-coded  mutation  probability 

0. 1  niching-parameter  for  constraints 


8.1.  Spatial  Discretization 

We  use  the  same  spatial  discretization  for  both  formulations.  For  the  confi  ned 
aquifer  we  discretize  the  domain  ft  =  [0, 1000]  x  [0, 1000]  x  [0,  30]  (m)  on  an 
equally  spaced  50  x  50  x  10  grid.  For  the  unconfi  ned  aquifer,  we  used  MOD- 
FLOW  to  determine  the  saturated  domain  ftunc  =  [0, 1000]  x  [0, 1000]  x 
[0, 27]  C  ft  and  then  discretized  ftunc  on  an  equally  spaced  50  x  50  x  10 
grid. 

8.2.  Five  well  formulation 
8.2.1.  Initial  Iterate 

IFFCO  requires  a  feasible  initial  iterate.  Figure  3  shows  the  steady  state  fbw 
fi  eld  for  the  confi  ned  aquifer.  Since  the  head  value  is  high  in  the  lower  left 
corner  we  initially  placed  one  well  there.  After  the  wells  are  activated,  the 
constraint  on  the  drawdown  is  violated  if  the  wells  are  too  close  together.  We 
looked  at  several  different  initial  iterates  until  we  found  one  that  satisfi  ed  the 
drawdown  constraint  for  both  the  confi  ned  and  unconfi  ned  aquifer.  We  found 
that  placing  the  remaining  four  wells  close  to  the  specifi  ed  head  boundaries 
and  signifi  candy  apart  from  each  other  was  feasible  for  both  physical  do¬ 
mains.  Figure  4  shows  the  relative  location  of  the  wells  and  the  pressure  head 
fi  eld  for  the  confi  ned  aquifer  with  the  wells  pumping  at  the  initial  iterate.  Note 
the  same  initial  well  locations  were  used  for  both  aquifers. 
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Figure  3.  Steady  state  head,  confined  aquifer 


Figure  4.  Initial  Iterate 
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We  refer  to  the  confi  ned  aquifer  as  CON  and  unconfi  ned  aquifer  as  UNC. 
The  function  value  at  the  initial  iterate  was  $23,204  for  CON  and  $26,958 
for  UNC.  Table  IV  shows  the  minimum  cost  found  and  the  number  of  calls 
to  MODFLOW  for  both  optimizers  and  aquifers.  IFFCO  reduced  the  cost  by 
6%  for  the  confined  aquifer  and  11%  for  the  unconfi  ned  aquifer.  For  both 
aquifers,  the  minimum  cost  found  by  the  GA  after  10  generations  was  5% 
higher  than  the  cost  found  by  IFFCO-which  is  high  considering  the  decrease 
from  the  initial  iterate.  Table  V  shows  the  initial  x-y  coordinates  for  the  5 
wells  and  the  optimal  locations  for  each  aquifer.  The  well  locations  at  the 
optimal  point  lie  on  the  boundary  constraint,  (15).  This  is  physically  rea¬ 
sonable,  since  the  head  values  arc  higher  in  that  region  due  to  the  Dirichlet 
boundary  conditions.  The  GA’s  cost  was  higher  than  IFFCO’s  because  well 
5  in  the  confi  ned  case  and  well  1  in  the  unconfi  ned  case  arc  not  relocated 
close  enough  to  the  Dirichlet  boundary  conditions  where  the  head  values  arc 
higher. 

In  our  evaluation  of  performance,  we  count  only  the  expensive  calls  to  the 
simulator  as  opposed  to  cumulative  calls  to  the  function.  This  is  the  approach 
taken  in  [6],  To  see  how  this  is  a  more  realistic  way  to  measure  cost,  consider 
the  case  where  the  linear  constraint  (12)  is  violated.  One  can  detect  this  viola¬ 
tion,  and  return  a  failure  for  /,  without  calling  the  expensive  fbw  simulator. 
IFFCO  is  being  modifi  ed  to  make  it  easy  for  the  user  to  evaluate  cost  at  a 
fi  ner  granularity  than  this,  to  allow  for  the  use  of  multiple  simulators  within 
the  evaluation  of  the  objective  function  and  constraints.  While  this  is  a  simple 
change  in  a  serial  code,  correctly  counting  the  calls  to  the  various  simulators 
in  a  parallel  implementation  requires  considerable  care. 

Figure  5  is  a  plot  of  the  value  of  the  objective  function  against  the  cumu¬ 
lative  number  of  calls  to  the  simulator.  IFFCO  is  currently  being  modifi  ed  to 
allow  the  user  to  easily  count  calls  to  the  objective  function  and  calls  to  the 
individual  simulators  that  arc  used  to  compute  it.  We  set  a  function  evaluation 
budget  of  10,000  for  this  work  and  IFFCO  converged  to  an  optimal  point 
within  approximately  3%  of  the  budget,  terminating  the  optimization  based 
on  the  sequence  of  finite  difference  scales.  Figure  5  shows  that  after  only 
roughly  100  function  evaluations,  the  objective  function  does  not  decrease 
signifi  candy. 

Figures  6  and  7  show  the  head  contours  in  the  layers  containing  the  wells 
with  the  wells  at  the  optimal  locations. 

8.3.  Six  well  lormulation 

The  results  above  arc  based  on  the  heuristic  that  installing  the  minimum 
number  of  wells  (5)  that  meet  the  extraction  target  is  the  best  approach.  To 
test  this,  we  compared  the  fi  ve  well  confi  guration  with  all  wells  pumping  at 
the  maximum  extraction  rate  to  a  six  well  confi  guration  with  both  locations 
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Table  IV.  Cost:  5  Wells 


Optimizer 

Problem 

min  f 

MODFLOW  Calls 

IFFCO 

CON 

$21,830 

275 

GA 

CON 

$22,822 

330 

IFFCO 

UNC 

$23,930 

302 

GA 

UNC 

$25,164 

328 

Table  V.  Optimal  Locations 


lnit_Co 

IFFCO 

GA 

IFFCO 

GA 

(m) 

CON  (m) 

CON  (m) 

UNC  (m) 

UNC  (m) 

X(l)  350.0 

401.7 

655.1 

464.2 

600.0 

Y(  1)725.0 

800.0 

737.8 

800.0 

216.2 

X(2)  775.0 

800.0 

794.3 

800.0 

397.5 

Y(2)  775.0 

800.0 

782.3 

800.0 

774.4 

X(3)  675.0 

776.9 

755.6 

800.0 

796.5 

Y(3)  675.0 

481.1 

203.1 

445.4 

706.0 

X(4)  200.0 

138.2 

569.3 

138.2 

149.7 

Y(4)  200.0 

800.0 

798.3 

800.0 

771.9 

X(5)  725.0 

798.4 

303.8 

800.0 

799.3 

Y(5)  350.0 

168.9 

501.6 

144.8 

513.3 

and  pumping  rates  as  decision  variables.  We  included  the  installation  cost  (/c 
in  (11))  in  the  objective  function  for  these  runs.  If  the  six  well  problem  is 
initialized  with  all  wells  pumping  at  the  maximum  extraction  rate,  then  one 
well  is  removed  from  the  design  in  the  course  of  the  optimization  and  the 
minimum  function  value  is  within  0.2%  of  that  found  with  the  original  fi  ve 
well  confi  guration.  If  the  six  wells  arc  initialized  with 

Qi  =  Qt  i  =  1  •  •  •  6, 

which  is  a  feasible  and  sensible  initial  iterate,  then  a  suboptimal  point  is 
found.  All  wells  remain  pumping  close  to  the  initial  pumping  rates,  although 
the  locations  align  with  the  specili  ed  head  boundary  conditions. 

The  objective  function  for  the  six  well  problem  contained  a  large  instal¬ 
lation  cost,  fc,  per  well.  A  common  approach  for  such  conditions  is  to  use 
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Figure  5.  Decrease  in  Function  Values 


Figure  6.  Confined  Aquifer 


a  mixed-integer  formulation  [37,  44,  30,  35],  Specifi  cally,  given  the  fact  that 
fc  was  signifi  candy  larger  than  f‘  and  that  a  minimum  of  fi  ve  wells  were 
required  to  satisfy  the  extraction  target,  a  reasonable  way  to  recast  the  prob¬ 
lem  was  to  include  an  integer  variable  indicating  which,  if  any,  well  should  be 
removed  from  the  design.  It  was  straightforward  to  include  an  integer  variable 
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s  ranging  from  1  to  8  in  the  GA’s  formulation.  A  value  of  s  in  the  range 
1, . . . ,  6  resulted  in  shutting  off  the  corresponding  well,  while  s  =  7, 8  led 
to  six  well  designs.  In  addition,  a  well  was  removed  from  the  design  if  its 
pumping  rate  fell  below  the  installation  threshold,  regardless  of  the  value 
of  s.  The  unequal  ranges  associated  with  fi  ve  and  six  well  designs  skewed 
the  GA’s  formulation  to  favor  fi  ve-well  designs.  This  refected  our  heuristic 
that  installing  the  minimum  number  of  wells  was  likely  to  be  cheaper  than  a 
six-well  design. 

Table  VI  shows  the  minimum  cost  and  the  number  of  calls  to  MODFLOW 
for  both  optimizers,  aquifers  for  the  better  initial  iterate.  The  cost  at  the  initial 
iterate  was  $170,972  for  the  confi  ned  aquifer  and  $152,878  for  the  unconfi  ned 
aquifer.  Both  the  GA  and  IFFCO  were  able  to  remove  one  well  from  the 
design,  resulting  in  a  solution  comparable  to  the  fi  ve  well  confi  guration  and 
reducing  the  cost  by  roughly  20%.  For  the  suboptimal  initial  iterate,  IFFCO 
was  only  able  to  decrease  the  cost  1%.  The  GA  did  not  find  anything  bet¬ 
ter  than  the  initial  iterate  in  30  generations  (over  100  function  evaluations) 
for  either  aquifer  with  the  suboptimal  initial  iterate  included  in  the  initial 
population. 

8.4.  Optimization  Landscapes 

The  fi  ve  well  confi  guration  does  not  have  a  discontinuous  installation  cost. 
To  get  a  better  understanding  of  the  objective  function  for  this  formulation, 
we  fi  xed  wells  2-5  and  computed  the  cost  while  letting  the  x  and  y  coordi- 
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Table  VI.  Cost:  6  Wells 


Optimizer 

Problem 

min  f 

MODFLOW  Calls 

1FFCO 

CON 

$140,237 

346 

GA 

CON 

$140,628 

464 

1FFCO 

UNC 

$124,582 

327 

GA 

UNC 

$127,069 

161 

nates  for  well  1  vary  between  20  and  800  meters.  Figure  8  shows  the  cost 
landscape  near  the  initial  iterate  for  the  confi  ned  aquifer  and  Figure  9  shows 
the  landscape  near  the  initial  iterate  for  the  unconfi  ned  aquifer.  The  peaks 
in  the  landscapes  occur  when  two  wells  get  close  together,  making  the  head 
values  low  and  hence  the  operational  cost  higher.  When  two  wells  get  too 
close  they  violate  (14),  leaving  a  small  infeasible  region  inside  each  of  the 
peaks.  These  peaks  also  make  the  landscapes  nonconvex  and  introduce  local 
minima.  When  we  tty  to  evaluate  the  function  at  an  infeasible  point,  we  do  not 
plot  an  ai'tifi  cial  value.  Note  that  only  a  subset  of  O  is  feasible,  especially  for 
the  unconfi  ned  aquifer  near  the  initial  iterate.  The  high  infeasibility  was  due 
to  repeated  violation  of  the  head  constraint  (14),  which  is  why  the  unconfi  ned 
case  is  more  challenging.  There  arc  also  small  discontinuities  apparent  in 
the  landscapes  since  we  round  real  numbers  to  grid  locations  to  run  the  fbw 
simulator. 

Figures  10  and  11  arc  the  surfaces  obtained  when  wells  2-5  arc  set  at  the 
optimal  locations  for  the  confi  ned  and  unconfi  ned  aquifers. 
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Figure  8.  Landscape  near  initial  iterate:  confined  aquifer 


xIO4 


Figure  9.  Landscape  near  initial  iterate:  unconfined  aquifer 
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Figure  10.  Landscape  near  solution:  confined  aquifer 
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Figure  11.  Landscape  near  solution:  unconfined  aquifer 
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8.5.  Discussion 

•  The  numerical  results  indicate  that  the  six  well  formulation  is  much  more 
diffi  cult  than  the  fi  ve  well  formulation  since  one  well  must  be  removed  from 
the  design  while  the  remaining  wells  extract  at  the  maximum  extraction  rate. 
IFFCO  did  well  on  the  the  six  well  formulation  with  an  initial  iterate  with  all 
wells  extracting  at  the  maximum  rate.  In  the  other  case  a  local  minimum  is 
found.  With  or  without  the  mixed-integer  formulation,  the  six  well  problem 
proved  diffi  cult  for  the  GA,  since  removing  a  well  only  led  to  a  feasible  iterate 
if  the  other  fi  ve  were  pumping  at  the  maximum  rate.  The  GA  was  able  to  fi  nd 
a  feasible  fi  ve  well  solution  when  the  better  initial  iterate  was  included  as  a 
member  of  the  initial  population.  With  a  completely  random  initial  population 
or  with  the  sub-optimal  initial  iterate  included,  the  GA  solution  for  the  con- 
fi  ned  aquifer  problem  was  a  sub-optimal  six  well  design.  The  same  was  true 
for  the  unconfi  ned  case.  However,  the  GA  with  a  completely  random  popu¬ 
lation  was  unable  to  fi  nd  a  feasible  iterate  for  the  unconfi  ned  problem  even 
when  using  a  population  size  of  two  hundred  and  running  for  two  hundred 
generations. 

•  We  also  considered  the  possibility  that,  over  a  longer  time  period,  the  six 
well  model  with  the  suboptimal  initial  iterate  may  be  superior  to  the  fi  ve  well 
model.  We  ran  both  the  confi  ned  and  unconfi  ned  problems  for  one  year  to 
determine  the  annual  operational  cost.  One  can  see  that  a  time  of  roughly 
130  years  for  the  confined  aquifer  and  90  years  for  the  unconfined  aquifer 
is  needed  to  obtain  a  lower  using  the  six  well  model.  Hence  the  fi  ve  well 
results  arc  the  most  realistic.  We  ran  both  problems  again  with  the  longer 
time  horizons  to  confi  rm  that  the  six  well  model  would  outperform  the  fi  ve 
well  model. 

•  It  is  common  in  practice  to  start  with  a  large,  fi  xed  grid  of  wells  and  seek 
only  the  optimal  pumping  rates,  removing  wells  from  the  design  as  needed. 
We  tried  this  approach,  using  N  =  16  so  that  we  seek  the  optimal  rates, 
{Qi}}=i  to  minimize  f7  =  fc  +  f°.  Neither  IFFCO  nor  the  GA  was  able 
to  converge  to  the  five  well  solution.  We  tried  several  different  initial  iter¬ 
ates  for  IFFCO  and  at  best,  in  429  function  evaluations,  IFFCO  had  left  10 
wells  in  the  design  and  reduced  the  cost  from  $  424,123  to  $  262,558  for 
the  confi  ned  aquifer.  The  results  were  similar  for  the  unconfi  ned  aquifer,  but 
it  was  even  more  challenging  to  fi  nd  a  feasible  initial  iterate.  We  performed 
several  experiments  with  IFFCO  and  with  a  good  initial  iterate,  IFFCO  was 
able  to  converge  to  the  fi  ve  well  solution  with  up  to  N  =  9  candidate  wells 
initially  fi  xed  on  the  grid.  With  a  relatively  large  number  of  fi  xed  wells,  it  was 
again  natural  to  use  a  mixed-integer  formulation  for  the  GA.  This  time,  the 
design  variables  consisted  of  a  rate  and  binary  “on-off”  switch  for  each  of  the 
A’  =  16  wells.  Starting  from  a  random  initial  population,  the  GA  produced  a 
six  well  design  for  both  the  confi  ned  and  unconfi  ned  problems.  The  fi  nal  cost 
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for  the  confined  aquifer  design  was  f  =  $161,588  after  510  function  eval¬ 
uations,  while  the  unconfi  ned  aquifer  solution  had  a  cost  of  f1  =  $142,755 
after  359  function  evaluations. 

•  This  application  is  challenging  for  formulations  with  N  >  5  wells  since 
wells  must  be  removed  from  the  design  space  to  decrease  the  installation  cost. 
Removing  wells  from  the  design  is  an  active  area  of  research  and  numerous 
approaches  exist  for  approximating  fixed  costs  with  continuous  functions 
(see  [32]  and  the  references  therein).  Our  approach  for  removing  a  well  if 
\Qi\  <  10  4  results  in  a  discontinuous  fi  xed  cost.  This  approach  was  imple¬ 
mented  due  to  its  simplicity,  but  also  because  the  implicit  fi  ltering  algorithm  is 
designed  for  problems  with  discontinuous  landscapes.  As  the  number  of  wells 
increased,  it  was  more  diffi  cult  to  fi  nd  a  feasible  initial  iterate  and  even  for  the 
six  well  formulation,  we  found  that  IFFCO  required  a  decent  initial  iterate.  A 
possible  remedy  for  this  might  be  to  use  another  optimization  routine  as  an 
initial  iterate  generator  for  IFFCO.  This  will  be  the  subject  of  future  work. 
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9.  Conclusions 

This  work  was  an  initial  analysis  of  a  subset  of  the  community  problems  pro¬ 
posed  in  [30].  The  formulation  of  the  objective  function  and  constraints  arc 
discontinuous  and  have  local  minima  as  the  optimization  landscapes  verify. 
These  features  arc  the  reason  for  the  failure  of  the  gradient-based  method 
(see  §  8).  As  pointed  out  in  [30],  deterministic  sampling  methods  have  not 
been  used  to  their  full  potential  in  the  subsurface  optimization  community. 
In  this  work  we  found  a  solution  to  the  well-fi  eld  design  application  with  the 
implicit  fi  ltering  algorithm  and  we  compared  our  results  to  those  obtained 
with  a  simple  genetic  algorithm. 

We  found  these  problems  to  be  challenging  for  several  reasons.  A  minimal 
cost  is  obtained  with  few  wells  pumping  at  large  rates.  For  these  problems, 
local  minima  exist  when  more  wells  arc  extracting  at  low  costs.  A  large 
decrease  in  the  objective  function  occurs  when  a  well  is  removed  from  the 
design  space.  The  installation  cost  for  this  work  is  discontinuous,  yet  we 
found  that  with  good  initial  data,  that  the  implicit  fi  ltering  algorithm  could 
perform  well  despite  the  discontinuous  formulation  for  N  <  9.  Another 
challenge  is  that  the  feasible  region,  especially  for  the  unconfi  ned  aquifer,  is 
small.  Although  implicit  fi  ltering  started  with  a  feasible  initial  iterate,  much 
experimental  work  was  done  to  fi  nd  one.  The  genetic  algorithm  was  unable 
to  fi  nd  a  feasible  point  for  the  six  well  formulation  and  required  a  feasible 
point  in  the  initial  population  in  order  for  the  optimization  to  progress  for  the 
unconfi  ned  aquifer. 

We  can  extend  this  study  to  improve  the  solution  for  this  type  of  problem. 

•  IFFCO  requires  a  feasible  initial  iterate,  and  the  numerical  results  show  that 
a  good  initial  iterate  is  needed  for  the  optimization.  Surrogate  models  based 
on  statistical  sampling  [5,  6]  may  be  a  good  way  to  explore  design  space  for 
good  initial  iterates. 

•  As  pointed  out  in  [30],  a  more  accurate  realization  of  the  subsurface  is 
needed  for  solutions  of  this  class  of  problems  to  be  used  in  decision  making. 
Adding  heterogeneities  to  the  domain  would  create  a  more  realistic  snap¬ 
shot  of  the  subsurface  yet  would  make  the  optimization  landscapes  much 
more  challenging.  More  robust  optimization  techniques  may  be  needed  as 
the  conceptual  domain  becomes  more  realistic. 

•  We  used  MODFLOW  to  simulate  fbw  for  a  well-fi  eld  design  application, 
despite  the  simulator’s  simple  well  model.  The  real-valued  well  location  that 
is  output  from  the  optimizer  is  rounded  to  a  grid  location.  For  a  more  accurate 
solution,  a  simulator  that  is  able  to  more  accurately  resolve  fbw  around  the 
well  is  essential.  A  well  model  that  need  not  place  wells  at  the  center  of  a  cell 
would  be  ideal. 

•  This  was  the  fi  rst  attempt  at  obtaining  a  solution  to  any  of  the  problems  pro¬ 
posed  in  [30].  An  in-depth  comparison  of  sampling  methods,  including  those 
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that  use  a  surrogate  response  surface,  is  currently  being  done  and  is  necessary 
before  any  solid  conclusions  can  be  made  on  which  method  performs  best  for 
this  class  of  problems. 


10.  Downloading  and  Running  the  Test  Problems 

The  problems  can  be  obtained  from 

http : //www4 . nesu . edu/ ~ctk/community . html 
The  test  problems  arc  packaged  as  compressed  UNIX  tar-  fi  les.  The  serial 
codes  arc  for  the  g77  compiler  and  have  been  tested  on  SUN  SparcStations 
running  Solaris,  various  Intel  platforms  running  Red  Hat  Linux  7.3  and  8.0, 
and  an  Apple  Macintosh  G4  running  OSX  10.2.  The  MPI  version  of  the  codes 
has  been  tested  on  an  IBM-SP3  and  a  DELL  Linux  server.  IFFCO  is  included 
in  the  packages.  The  README  fi  les  in  the  main  directory  explain  how  to 
assemble  the  fi  les  and  interpret  the  results. 

MODFLOW  can  be  obtained  directly  from  the  USGS  at  the  URL 
http: //water.usgs. gov/software/modf low- 9  6 . html 
The  USGS  provides  compiled  executables  for  SUN,  SGI,  and  DOS  sys¬ 
tems,  as  well  as  UNIX  source.  Our  packages  provide  makefi  les  for  some  other 
UNIX  environments. 
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